Tracking the Cracking: A Holistic Analysis of Rapid Ice Shelf Fracture Using Seismology, Geodesy, and Satellite Imagery on the Pine Island Glacier Ice Shelf, West Antarctica

Abstract Ice shelves regulate the stability of marine ice sheets. We track fractures on Pine Island Glacier, a quickly accelerating glacier in West Antarctica that contributes more to sea level rise than any other glacier. Using an on‐ice seismic network deployed from 2012 to 2014, we catalog icequakes that dominantly consist of flexural gravity waves. Icequakes occur near the rift tip and in two distinct areas of the shear margin, and TerraSAR‐X imagery shows significant fracture in each source region. Rift‐tip icequakes increase with ice speed, linking rift fracture to glaciological stresses and/or localized thinning. Using a simple flexural gravity wave model, we deconvolve wave propagation effects to estimate icequake source durations of 19.5–50.0 s and transient loads of 3.8–14.0 kPa corresponding to 4.3–15.9 m of crevasse growth per icequake. These long‐source durations suggest that water flow may limit the rate of crevasse opening.

X -2 : dynamics. After obtaining the 30 s ITRF solutions, we performed a 5 min weighted average using the inverse of the individual epoch uncertainties for data weights, and then rotated the XYZ displacements into local North, East, and up displacements.
We obtain ice speed from the processed GPS positions at the GPS station SOW3 by calculating the total distance moved in each day of the deployment and differentiating with respect to time. The resulting ice speed curve contains some spike artifacts that arise from numerical differentiation, which we remove by linearly interpolating between the ice speed before and after the affected time period. Finally, we low pass filter the data to remove trends on time periods shorter than a week.

Icequake detection
To detect flexural gravity icequakes in the dataset, we design a two-stage detection scheme that identifies broadband, dispersive seismic events. First, we employ a short term average/long term average (STA/LTA) impulsivity detector. This method identifies high-amplitude impulsive events by comparing the mean amplitude of a short time window with the mean amplitude of a long time window (Allen, 1978). The detector is triggered when STA exceeds LTA by some threshold. STA/LTA threshold values are selected by tuning the algorithm to successfully detect high signal-to-noise ratio manually-identified events (see Table S1). We carry out STA/LTA on the vertical component of each station separately in two different frequency bands (0.01-1 Hz and 1-10 Hz). Selected waveforms satisfy the STA/LTA trigger criteria in both frequency bands on at least three out of the five stations. We refine the catalog and generate waveform templates by cross-correlating April 15, 2022, 7:09pm : X -3 each preliminary event with a master event waveform and selecting the events with cross correlation coefficients exceeding 0.9. This selection procedure resulted in 57 template events.
Second, we perform a template matching technique based on cross-correlation to identify events that were similar to the events in the preliminary catalog (Gibbons & Ringdal, 2006). To detect new events, each template event is cross correlated with all time windows in the dataset on two frequency bands (0.05-1 Hz and 1-10 Hz). We increase the lower frequency bound from 0.01 Hz to 0.05 Hz since many template events contained uninterpretable noise at frequencies below 0.05 Hz. The detector is triggered when the cross-correlation coefficient between a template event waveform and the given time window exceeds a threshold. The threshold value is selected so that the algorithm successfully detects the other known events of the preliminary catalog (see Supporting Table S1). Detected waveforms satisfy the trigger criteria on at least three out of the five stations in both frequency bands. We carry out this procedure for each template and removed redundant detections to yield the final catalog.
We detect 22,119 seismic events using the two-band template matching scheme. The detected events have a typical duration of around 50 s and an average peak vertical velocity of approximately 1 · 10 −5 m/s. Event waveforms vary in shape, indicating varied sources and propagation paths. Many of the events exhibit characteristic dispersion between 0.05 and 1 Hz with high frequencies arriving before low frequencies, while others were monochromatic between 0.05 and 1 Hz.

Waveform Clustering
Because the catalog of detected events contains both dispersive and monochromatic waveforms, we seek to cluster the events into groups based on wave shape. To do so, we modify the K-shape algorithm of Paparrizos and Gravano (2016). K-shape is designed specifically to cluster time series data. Instead of calculating the Euclidean distance between potential cluster centers and observations, K-shape calculates distances using the maximum normalized cross correlation coefficient between two time series. We adapt the K-shape algorithm for three component seismic data by independently computing the cross-correlation time series between the three separate seismic channels (vertical, East, and North). We then sum these three cross-correlation time series and calculate the distance metric as the maximum value of this summed cross correlation time series.
We use the K-shape algorithm to divide the catalog into 2, 3, · · · , 20 clusters. However, beyond two clusters, the differences between waveforms in each cluster become progressively less clear, and an analysis of the average distance from waveforms to their cluster center does not show significant improvement for larger numbers of clusters. We thus use the K-shape algorithm to divide the catalog into two distinct clusters, which differ based on waveform dispersion. The first cluster contains 8,184 dispersive events. The second cluster contains 13,935 monochromatic events that do not exhibit dispersion within the chosen frequency band. This difference suggests that the two types of waveforms may have been generated by different source processes. Since we are specifically interested in dispersive flexural gravity wave signals, we restrict the remaining analysis to the dispersive cluster. : X -5 3. Text S3. Methods for computing event back-azimuths.

Robust first arrival determination
We obtain the relative first arrival time of each event through phase lags measurements.
We cross-correlate each respective component waveform between each seismic station.
We choose a window length of 500 s around the first arrival. The trace that requires the largest shift forward in time to align with the other traces is taken to be the station of first arrival. In most cases, the first arrivals obtained independently using each component are in agreement for at least two components out of three. However, if all three components produce different stations of first arrival, a back-azimuth is not calculated and the event is disregarded.

Amplitude threshold
Next, we ensure that the polarization is extracted over a high signal to noise ratio event as against noise. We slide through the event waveform in 10 s windows with a step size of five seconds. For each 10 s time window, we check if the average amplitude of that window exceeds the average amplitude of the entire 500 s event window.

Principal component analysis
For time windows with sufficiently large amplitude, principal component analysis (PCA) is performed on the HHE (East) and HHN (North) traces from each station to retrieve the PCA components. The PCA first component is a vector whose direction explains the largest contribution of the data variance. It is equivalent to the eigenvector of the data covariance matrix that has the largest eigenvalue.

Determining the predicted first arrival
We try three methods of computing the predicted station of first arrival corresponding to both possible propagation directions.
In the first method, we compare both possible phase back-azimuths to the back-azimuths of each station with respect to the mean station location, or array centroid. The stations that are radially closest to each possible back-azimuth are predicted to be the two possible first arrivals. The sign of the PCA first component vector is then adjusted to match the propagation direction whose predicted first arrival agree with the observed first arrival.
Phases for which neither predicted first arrival agreed with the observed first arrival are discarded.
In the second method, we divide the array into two sectors along a line through the array centroid orthogonal to the PCA first component vector. The sign of the PCA first April 15, 2022, 7:09pm : X -7 component vector is then adjusted to match the propagation direction corresponding to the sector containing the observed first arrival. No phases are discarded.
In the third method, we compute the distance vector from the array centroid to each station. For incoming plane waves, the station farthest from the array centroid in the direction of propagation records the first arrival. The stations whose distance vectors have the largest component oriented in each possible propagation directions are predicted to be the two possible first arrivals. The sign of the PCA first component vector is then adjusted to match the propagation direction whose predicted first arrival agree with the observed first arrival. Phases for which neither predicted first arrival agreed with the observed first arrival are discarded. All three methods gave relatively consistent results.

Back-azimuth stacking
Next, we sum the PCA first component vectors across each station to obtain an average vector whose norm indicates the level of agreement between propagation directions calculated at each station. Finally, we take the arctangent of the quotient of the two elements of the PCA component vector to retrieve a back-azimuth. Because this procedure is repeated for each 10 s time window in the event, the result for each individual event is a distribution of back-azimuths calculated for each time window within that event.
To obtain a single back-azimuth for each event, we take the average of the back-azimuths calculated using each time window in the data. We use the mean of circular quantities, with the back-azimuth from each time window weighted by the norm of the summed PCA components across the array for that window. This means that time windows with poor agreement between stations are downweighted when taking the average back-azimuth. The April 15, 2022, 7:09pm X -8 : weighted mean of circular quantities is expressed below for the back-azimuth distribution θ 1 , ..., θ n with PCA norms w 1 , ..., w n of a single event with n time windows:

Uncertainties in icequake back-azimuths
Uncertainties for icequake back-azimuths were computed using the weighted standard deviation of circular quantities for each event. The weighted standard deviation of circular quantities is expressed below for the back-azimuth distribution θ 1 , ..., θ n with PCA norms w 1 , ..., w n of a single event with n time windows: The mean back-azimuthal uncertainty for rift-tip icequakes was 26.16 suggests that the recorded icequakes were generated by that fracture. Third, seismicity in front of a rift's visible extent was observed at the Amery Ice Shelf by Bassis et al. (2007) and may indicate that the wing crack was more advanced at depth than at the surface, explaining events whose back-azimuths point just in front of the crack. Finally, no other areas consistent with distribution of rift-tip icequake back-azimuths contain any visible fracture opening or growth, suggesting that the large amount of fracture evolution at the rift tip observed in imagery is the most likely source of the recorded icequakes.

Analytical Solution for Ocean Surface Waves
We examine the water velocity potential function ϕ and relate it to the vertical ice shelf velocity w. We first solve the ocean surface wave equation for a body of water with infinite length and finite depth: over the interval −∞ < x < ∞, −h w < y < 0. We enforce zero velocity at the ocean floor and couple vertical velocity to the rate of beam deflection at the ocean surface: We enforce the Sommerfeld radiation condition: We employ the Fourier Transform, written for an arbitrary function f (x) as Applying the Fourier Transform in x to each term of Equation 3 yields: Evaluating the integral term: We then apply the radiation condition (Equation 5) to eliminate the first two terms, resulting in an ordinary differential equation ofφ: By applying the vertical boundary conditions (Equation 4), we obtain the timewavenumber domain solution that satisfies the governing equation and boundary con- We note that ϕ is a linear function of w, therefore permitting us to write the floating beam equation using the linear operator A as noted in the main text.

Analytic Solution for Buoyant Ice Shelf Flexure
To interrogate the source process that explains the observations, we obtain the Green's function, or fundamental solution of a floating dynamic beam to an impulse forcing. We April 15, 2022, 7:09pm : X -11 obtained the Green's function by using integral transform methods to solve the governing equation for an impulse forcing in space and time. We write the Green's function formulation of Equation 1 from the main text: where G is the Green's function, δ(x) is Dirac delta function in space, and δ(t) is the Dirac delta function in time. As before, we apply the Fourier Transform in space to each term. Next, we apply the Laplace transform, defined as, Applying both the Fourier and Laplace Transforms to Equation 8 yields: We can then solve forḠ * algebraically: Finally, we analytically compute the inverse Laplace transform of Equation 9 to obtain the Fourier-transformed Green's function, In practice, we numerically calculateḠ for a range of times and wavenumbers that define the temporal and spatial domain of the model run. OnceḠ is calculated for each element of a vector of times and a vector of wavenumbers, the IFFT (inverse fast Fourier X -12 : transform) is taken to numerically retrieve the Green's function G(x, t) of the ice shelf for an applied unit point force.

Greens function for a point moment source
To retrieve the impulse response to a point bending moment source, we note that an applied bending moment is equivalent to a pair of infinitesimally-spaced point loads with opposite signs: To obtain G m (x, t), we numerically take the spatial derivative of the point load Green's function G(x, t).

Deconvolution procedure
We calculate source load through the deconvolution, where hats denote Fourier-transformed quantities, F −1 is the inverse Fourier transform, w observed (t) is a linear stack of observed displacement seismograms, P estimated (t) is an estimated source load distribution, and x 0 is the station epicentral distance. We obtain w observed (t) for each spatial group by aligning each waveform in the group with respect to a master event using cross correlation and taking the average waveform. Master events were selecting by finding the event from each spatial group that was best-correlated with April 15, 2022, 7:09pm : X -13 the overall centroid of the dispersive cluster. We choose x 0 corresponding to the average distance to each spatial group: for the rift tip, x 0 = 25 km; for rift/margin, x 0 = 25 km; for margin icequakes, x 0 = 17.5 km. We alternatively consider a bending moment source through the relationship, 5. Text S5, Estimating vertical crack extent from point load magnitude The point load applied during vertical crack growth arises from the difference in ice overburden stress above the crack and buoyancy stress exerted by the water that fills the crack. Therefore, we write the following expression and substitute in the each stress: where P is the applied point load, ρ i is the density of ice, ρ w is the density of water, h i is the ice shelf thickness, ∆z c is the change in vertical position of the crack tip, and z w is the position of the water line.
Substituting the expression for the water line position z w = ρ i ρw h i and simplifying gives: Finally, we write: Table S1.

References
April 15, 2022, 7:09pm   have larger amplitude and longer duration for thicker beams, because larger forcing is required to induce a given displacement for a more rigid beam. Flexural rigidity, the parameter that governs flexure, is a function of thickness.
Modeled beam thicknesses are shown in the legend. Source time functions generally have larger amplitude and longer duration for thicker beams, because larger forcing is required to induce a given displacement for a more rigid beam. Flexural rigidity, the parameter that governs flexure, is a function of thickness.
April 15, 2022, 7:09pm  The resulting modeled displacements, shown in black, have a longer decay and larger amplitdue low-frequency displacements than the rift/margin stack, shown in green, for both bending moment and point load sources.
The resulting modeled displacements, shown in black, have a longer decay and larger amplitdue low-frequency displacements than the shear margin stack, shown in purple, for both bending moment and point load sources. The modeled displacements arising from an applied bending moment are relatively similar to the shear margin stack, but the results of deconvolution do not support the hypothesis that the observations were generated by a step forcing in bending moment.